rm(list=ls())

# set WD here

library(stargazer)
library(ggplot2)
library(dplyr)
library(lfe)
library(psych)


## Read in data:
first_experiment <- read.csv("first_experiment.csv") #data from the first experiment (2018)
first_experiment$exp <- 1

second_experiment <- read.csv("second_experiment.csv") #data from the second experiment (2019)
second_experiment$exp <- 2


## Merge the data:
combined <- rbind(first_experiment, second_experiment)


## Some data cleaning:
combined$station1 <-
  gsub(" ", "", combined$station1, fixed = TRUE) #remove empty spaces from station names


## Creating some additional variables (binary variables for heat (at different thresholds)):
combined$hot24 <- ifelse(combined$temp>24,1,0)
combined$hot25 <- ifelse(combined$temp>25,1,0)
combined$hot26 <- ifelse(combined$temp>26,1,0)

combined_hot24 <- combined
combined_hot24$temperature <- combined_hot24$hot24

combined_hot25 <- combined
combined_hot25$temperature <- combined_hot25$hot25

combined_hot26 <- combined
combined_hot26$temperature <- combined_hot26$hot26



####################################
# TABLES
####################################

############
# Table 1
# Balance Tests: Bystander Characteristics: Hot vs. Normal Temperature
############

balance <- 
  rbind(
    c("Share of women bystanders",
      summary(felm(sharewomen~temperature | rush + station1 + bystander  
                   | 0 | station1, data=combined_hot25))$coefficients[1,][c(1,2,4)]),
    c("Share of bystanders with earphones",
      summary(felm(shareearphones~temperature | rush + station1 + bystander  
                   | 0 | station1, data=combined_hot25))$coefficients[1,][c(1,2,4)]),
    c("Share natives bystanders",
      summary(felm(sharenative~temperature | rush + station1 + bystander  
                   | 0 | station1, data=combined_hot25))$coefficients[1,][c(1,2,4)]),
    c("Share of bystanders below 30",
      summary(felm(sharebelow30~temperature | rush + station1 + bystander  
                   | 0 | station1, data=combined_hot25))$coefficients[1,][c(1,2,4)]),
    c("Share of bystanders above 60",
      summary(felm(shareabove60~temperature | rush + station1 + bystander  
                   | 0 | station1, data=combined_hot25))$coefficients[1,][c(1,2,4)]),
    c("Share of Christian bystanders",
      summary(felm(christian~temperature | rush + station1 + bystander  
                   | 0 | station1, data=combined_hot25))$coefficients[1,][c(1,2,4)]),
    c("Share of non-religious bystanders",
      summary(felm(noreligion~temperature | rush + station1 + bystander  
                   | 0 | station1, data=combined_hot25))$coefficients[1,][c(1,2,4)]),
    c("Share of bystanders full-time employed",
      summary(felm(workft~temperature | rush + station1 + bystander  
                   | 0 | station1, data=combined_hot25))$coefficients[1,][c(1,2,4)]),
    c("Share of bystanders with university education",
      summary(felm(uniplus~temperature | rush + station1 + bystander  
                   | 0 | station1, data=combined_hot25))$coefficients[1,][c(1,2,4)])
  )

colnames(balance) <- c("Variable","Difference","SE", "p-value")
balance[,2] <- round(as.numeric(balance[,2]),digits=3)
balance[,3] <- round(as.numeric(balance[,3]),digits=3)
balance[,4] <- round(as.numeric(balance[,4]),digits=3)

stargazer(balance, title="Balance Tests", type = "text")
#stargazer(balance, title="Balance Tests")



############
# Table 2
# Help Behavior by Temperature
############

temperature_absolute <- felm(anyhelp~temp*treat | 0 | 0 | station1, data=combined)
temperature_absolute_fe <- felm(anyhelp~temp*treat | rush + station1 | 0 | station1, data=combined)
temperature_absolute_fe_bystander <- felm(anyhelp~temp*treat | rush + station1 + bystander | 0 | station1, data=combined)

stargazer(temperature_absolute, temperature_absolute_fe, temperature_absolute_fe_bystander, 
          title="Help Behavior by Temperature", align=TRUE,star.cutoffs = c(0.2, 0.1, 0.02), type = "text", notes = "*p<0.1; **p<0.05; ***p<0.01, one-tailed test.",notes.append = FALSE)
#stargazer(temperature_absolute, temperature_absolute_fe, temperature_absolute_fe_bystander, 
#          title="Help Behavior by Temperature", align=TRUE,star.cutoffs = c(0.2, 0.1, 0.02), notes = "Estimated with linear regression. Standard errors (clustered at the station level) in
#          parentheses. *p<0.1; **p<0.05; ***p<0.01, one-tailed test.",notes.append = FALSE)



############
# Additonal Tests Referenced in the Manuscript
############

# interaction between temperature, discrimination, and rush hour
# to probe whether there is heterogeneity in the relationship between temperature and discrimination 
# between iterations conducted during rush hours and other times of the day
# (discussed on page 9):
summary(felm(anyhelp~temp*treat*rush | station1 | 0 | station1, data=combined))


# help behavior by temperature toward outgroup members only
# (discussed on page 13):
summary(felm(anyhelp~temp | 0 | 0 | station1, data=subset(combined, treat==1)))



############
# Table S1
# Descriptive Statistics
############

combined_valid <- subset(combined,!is.na(anyhelp) & !is.na(temp) & (treat==1 | treat==0)) #subsetting to valid observations on key variables
descript <- combined_valid %>% 
  select(temp, anyhelp, sharewomen, shareearphones, sharenative, 
         sharebelow30, shareabove60, christian, noreligion, workft, uniplus)

descriptives <- describe(descript)[,c(8,9,3,5,4,2)]
descriptives <- 
  cbind(c("Temperature","Assistance (Outcome)","Share of women bystanders",
          "Share of bystanders with earphones", "Share native bystanders", "Share of bystanders below 30", 
          "Share of bystanders above 60", "Share of Christian bystanders", "Share of non-religious bystanders", 
          "Share of bystanders full-time employed", "Share of bystanders with university education"),descriptives)
names(descriptives) <- c("Variable","Min","Max","Mean","Median","SD","N")
descriptives$Mean <- round(descriptives$Mean,2)
descriptives$SD <- round(descriptives$SD,2)

stargazer(as.matrix(descriptives),rownames=F, title="Descriptive Statistics", type="text")
#stargazer(as.matrix(descriptives),rownames=F, title="Descriptive Statistics")



############
# Table S2
# Help Behavior at Hot Temperatures
############

results_hot24 <- felm(anyhelp~temperature*treat | 0 | 0 | station1, data=combined_hot24)
results_hot24_fe <- felm(anyhelp~temperature*treat | rush + station1 | 0 | station1, data=combined_hot24)

results_hot25 <- felm(anyhelp~temperature*treat | 0 | 0 | station1, data=combined_hot25)
results_hot25_fe <- felm(anyhelp~temperature*treat | rush + station1 | 0 | station1, data=combined_hot25)

results_hot26 <- felm(anyhelp~temperature*treat | 0 | 0 | station1, data=combined_hot26)
results_hot26_fe <- felm(anyhelp~temperature*treat | rush + station1 | 0 | station1, data=combined_hot26)

stargazer(results_hot24, results_hot24_fe, results_hot25, results_hot25_fe, results_hot26, results_hot26_fe, type="text", title="Help Behavior at Hot Temperatures", align=TRUE,star.cutoffs = c(0.2, 0.1, 0.02), notes = "*p<0.1; **p<0.05; ***p<0.01, one-tailed test.",notes.append = FALSE)

#stargazer(results_hot24, results_hot24_fe, results_hot25, results_hot25_fe, results_hot26, results_hot26_fe, title="Help Behavior at Hot Temperatures", align=TRUE,star.cutoffs = c(0.2, 0.1, 0.02), notes = "Estimated with linear regression. Standard errors (clustered at the station level) in
#parentheses. *p<0.1; **p<0.05; ***p<0.01, one-tailed test.",notes.append = FALSE)



############
# Table S3
# Help Behavior by Temperature: East vs. West
############

temperature_absolute_region <- felm(anyhelp~temp*treat*east | 0 | 0 | station1, data=combined)
temperature_absolute_fe_region <- felm(anyhelp~temp*treat*east | rush | 0 | station1, data=combined)
temperature_absolute_fe_bystander_region <- felm(anyhelp~temp*treat*east | rush + bystander | 0 | station1, data=combined)

stargazer(temperature_absolute_region, temperature_absolute_fe_region, temperature_absolute_fe_bystander_region, type = "text",
          title="Help Behavior by Temperature: East vs. West", align=TRUE,star.cutoffs = c(0.2, 0.1, 0.02), notes = "*p<0.1; **p<0.05; ***p<0.01, one-tailed test.",notes.append = FALSE)


#stargazer(temperature_absolute_region, temperature_absolute_fe_region, temperature_absolute_fe_bystander_region, 
#          title="Help Behavior by Temperature: East vs. West", align=TRUE,star.cutoffs = c(0.2, 0.1, 0.02), notes = "Estimated with linear regression. Standard errors (clustered at the station level) in
#          parentheses. *p<0.1; **p<0.05; ***p<0.01, one-tailed test.",notes.append = FALSE)



####################################
# FIGURES
####################################

############
# Figure 2
# Help rates to natives or immigrants with hijab by temperature level with linear trendlines
############

fig2 <- 
  ggplot(combined_valid, aes(x = temp, y = anyhelp, color=factor(treat))) + 
  geom_point(size = 3, alpha = 0.3) + theme_minimal() +
  stat_smooth(method = "lm", formula = y ~ poly(x, 1), se = TRUE, level=0.95, span = 0.9) + 
  xlim(17, 35) +
  ylim(0, 1) +
  xlab("Temperature (Celsius)") +
  ylab("Assistance Rates") +
  scale_colour_manual(name="Treatment",
                      values=c("#0072B2","#CC0000"),
                      labels = c("Native", "Immigrant\nwith Hijab")) 
suppressWarnings(ggsave("Figure2.pdf", fig2, width = 7, height = 6, units = "in"))
suppressWarnings(print(fig2, width = 7, height = 6, units = "in"))



############
# Figure S5
# Help rates in response to the two treatments by absolute temperature with LOESS curves
############

figS5 <- 
  ggplot(combined_valid, aes(x = temp, y = anyhelp, color=factor(treat))) + 
  geom_point(size = 3, alpha = 0.3) + theme_minimal() +
  stat_smooth(method = "loess", formula = y ~ poly(x, 1), se = TRUE, level=0.95, span = 0.9) + 
  xlim(17, 35) +
  ylim(0, 1) +
  xlab("Temperature (Celsius)") +
  ylab("Assistance Rates") +
  scale_colour_manual(name="Treatment",
                      values=c("#0072B2","#CC0000"),
                      labels = c("Native", "Immigrant\nwith Hijab")) 
suppressWarnings(ggsave("FigureS5.pdf", figS5, width = 7, height = 6, units = "in"))
suppressWarnings(print(figS5, width = 7, height = 6, units = "in"))

